/*
 * get_xplor.c
 *
 *  Created on: 21 avr. 2010
 *      Author: marchi
 */
#include "stdio.h"
#include "stdlib.h"
#include "math.h"
#include "get_xplor.h"
#include "physics.h"

/* WARNING:
 * Do _not_ use these routines to calculate the angle between two vectors
 * as acos(cos_angle(u,v)). While it might seem obvious, the acos function
 * is very flat close to -1 and 1, which will lead to accuracy-loss.
 * Instead, use the new gmx_angle() function directly.
 */
static inline real
cos_angle_no_table(const rvec a,const rvec b)
{
  /* This version does not need the invsqrt lookup table */
  real   cosval;
  int    m;
  double aa,bb,ip,ipa,ipb; /* For accuracy these must be double! */

  ip=ipa=ipb=0.0;
  for(m=0; (m<DIM); m++) {		/* 18		*/
    aa   = a[m];
    bb   = b[m];
    ip  += aa*bb;
    ipa += aa*aa;
    ipb += bb*bb;
  }
  cosval=ip/sqrt(ipa*ipb); 		/* 12		*/
					/* 30 TOTAL	*/
  if (cosval > 1.0)
    return  1.0;
  if (cosval <-1.0)
    return -1.0;

  return cosval;
}

static inline real norm2(const rvec a)
{
  return a[XX]*a[XX]+a[YY]*a[YY]+a[ZZ]*a[ZZ];
}
static inline real norm(const rvec a)
{
  return (real)sqrt(a[XX]*a[XX]+a[YY]*a[YY]+a[ZZ]*a[ZZ]);
}

static int cell_angles(matrix box, real *a, real *b, real *g){
	real alpha,beta,gamma;
	if (norm2(box[YY])*norm2(box[ZZ])!=0)
		alpha = RAD2DEG*acos(cos_angle_no_table(box[YY],box[ZZ]));
	else
		alpha = 90;
	if (norm2(box[XX])*norm2(box[ZZ])!=0)
		beta  = RAD2DEG*acos(cos_angle_no_table(box[XX],box[ZZ]));
	else
		beta  = 90;
	if (norm2(box[XX])*norm2(box[YY])!=0)
		gamma = RAD2DEG*acos(cos_angle_no_table(box[XX],box[YY]));
	else
		gamma = 90;
	*a=alpha;
	*b=beta;
	*g=gamma;
	return 0;
}

int get_xplor(real *** Phi,int nx, int ny, int nz, matrix box, FILE *fp){
	int nx_s,ny_s,nz_s;
	int nx_e,ny_e,nz_e;
	real a=0.0,b=0.0,c=0.0,alpha=0.0, beta=0, gamma=0;

	nx_s=0;
	ny_s=0;
	nz_s=0;
	nx_e=nx-1+nx_s;
	ny_e=ny-1+ny_s;
	nz_e=nz-1+nz_s;

	cell_angles(box, &alpha, &beta, &gamma);
	a=norm(box[XX]);
	b=norm(box[YY]);
	c=norm(box[ZZ]);

	fprintf(fp,"\n       2 !NTITLE\n") ;
	fprintf(fp," REMARKS Electrostatic Potential from GROMACS\n") ;
	fprintf(fp," REMARKS DATE: 2010-04-21\n") ;
	fprintf(fp," %7d %7d %7d %7d %7d %7d %7d %7d %7d\n",nx,nx_s,nx_e,ny,ny_s,ny_e,nz,nz_s,nz_e);
	fprintf(fp,"%12.5E%12.5E%12.5E%12.5E%12.5E%12.5E\n",10*a,10*b,10*c,alpha,beta,gamma);
	fprintf(fp, "ZYX\n") ;
	int ja;
	for(int k=0;k<nz;k++){
		fprintf(fp,"%8d\n",k);
		ja=0;
		for(int j=0;j<ny;j++){
			for(int i=0;i<nx;i++){
				ja++;
				fprintf(fp,"%12.5E",Phi[i][j][k]);
				if(!(ja%6)) fprintf(fp,"\n");
			}
		}
		fprintf(fp,"\n");
	}
	fprintf(fp,"%8d\n",-9999);
	return 0;
}
